Reconstructing the quantum state of oscillator networks with a single qubit 
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We introduce a scheme to reconstruct arbitrary states of networks composed of quantum 
oscillators — e.g., the motional state of trapped ions or the radiation state of coupled cavities. The 
scheme involves minimal resources and minimal access, in the sense that it i) requires only the 
interaction between a one-qubit probe and a single node of the network; ii) provides the Weyl char- 
acteristic function of the network directly from the data, avoiding any tomographic transformation; 
m) involves the tuning of only one coupling parameter. In addition, we show that a number of 
quantum properties can be extracted without full reconstruction of the state. The scheme can be 
used for probing quantum simulations of anharmonic many-body systems and quantum computa- 
tions with continuous variables. Experimental implementation with trapped ions is also discussed 
and shown to be within reach of current technology. 
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Coupled harmonic and anharmonic oscillators consti- 
tute the building blocks of mathematical models that are 
ubiquitous in physics. Quantum systems are no excep- 
tion and, in fact, the studies on coupled quantum os- 
cillators trace back to the origin of quantum physics it- 
self. The emergence of quantum information science has 
added renewed interest in these continuous- variable sys- 
tems For example, the possibility to exert exquisite 
experimental control over travelling oscillator modes led 
to novel applications in quantum optical communication 
0. However, these experiments involve solely a lim- 
ited number of modes, whereas both the investigation 
on many-body models and the quest for advanced quan- 
tum information tasks call for the realization of more 
complex bosonic networks. Interestingly, some alterna- 
tive experimental settings are now reaching maturity for 
implementing these networks, thanks to unprecedented 
ability to manipulate confined quantum modes Q. In 
fact, trapped ions, cavity QED, circuit QED, or nanome- 
chanical oscillators have been proposed to realize quan- 
tum simulators of many-body systems whose properties 
are beyond reach of purely theoretical and numerical in- 
vestigations Suggestions for the use of these plat- 
forms for continuous-variable quantum computation [sj 
have also been recently put forward @ . In addition, gen- 
eral physical concepts, related to entropy-area laws 
or quantum thermodynamics [1], have been extensively 
analyzed for oscillator networks and could be amenable 
for experimental testing. 

Despite the aforementioned proposals, fundamental 
tools that still lack in this context are minimal and fea- 
sible schemes to reconstruct the quantum state of oscil- 
lator networks — a necessary step for probing the validity 
of quantum simulations and computations. In general, 
a reconstruction scheme — also dubbed quantum state 
tomography — tries to estimate a quantum state using 
measurements on an ensemble of identical copies of it. 
Considering travelling modes a huge research effort has 
been made in the past years and quantum tomography 
is now standard @. However, the latter is based on the 
measure of quadrature signals, which are unavailable for 
confined quantum modes. To face this obstacle, in the 



case of a single oscillator, many alternative schemes have 
then been put forward, relying on interrogating the sys- 
tem either with a discrete- variable (qubit) jlOl [T]| or a 
continuous-variable [12| probe. However, the adaptation 
of these schemes to the relevant case of a network of many 
oscillators has been vastly overlooked (see Refs. [l^ for 
some details). In quantum tomography of travelling op- 
tical fields, the state reconstruction of a iV-modc field re- 
quires the ability to perform joint measurements of N ar- 
bitrary field quadratures (one for each mode) . Similarly, 
if we wanted to reconstruct the joint state of N oscilla- 
tors with a probe-mediated method, we might introduce 
N auxiliary probes. This approach, requiring maximal 
access to the network, can quickly become impractical 
as the number of oscillators is increased. In fact, ac- 
cessibility constraints often plague experiments and thus 
the question of extracting information with only partial 
access — a non-trivial inverse problem involving an inter- 
acting many-body system — assumes also a practical rele- 
vance. For example, two recent experiments with trapped 
ions have probed the dynamics of the simplest possible 
network, composed of two oscillators [iJl- There, only 
one of the oscillators could be probed, imposing partial 
access to the system. In general, it is thus desirable to de- 
sign state reconstruction protocols that involve a smaller 
number of resources, as compared to the straightforward 
extension of the single-oscillator schemes. We introduce 
here one such protocol, that solves these major draw- 
backs by requiring only minimal access to the network. 
In particular, it involves only the interaction between one 
qubit-probe and one constituent of the network. In addi- 
tion, the method provides directly the Weyl character- 
istic function of the system, avoiding the massive post- 
processing of noisy data common to many reconstruction 
schemes — a benefit that considerably eases its implemen- 
tation. 

We consider a generic oscillator network in an unknown 
state — possibly being an eigenstate of some simulated an- 
harmonic model, or an intermediate state of a quantum 
computation. Regardless the previous dynamics, we sup- 
pose that, from a certain time t = 0, the oscillators in- 
teract only harmonically. In addition, a single qubit can 
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interact with a single (fixed) oscillator via a tunable bi- 
linear coupling (see Fig[T]). Such a network-probe system 
is then let evolve for a certain period of time, allowing 
part of the information about the network state to be 
transfered into the qubit. Afterwards, only the qubit is 
measured. Repeating the procedure it is possible to re- 
construct the state of the whole network, by solely tuning 
the profile of the interaction strength. 

The paper is structured as follows. In Sec. |T] we intro- 
duce the Hamilonian model and solve the time dynamics 
of the system. We will see in Sec. |lT] that the system 
dynamics performs arbitrary qubit-controlled multimode 
displacements of the network. This, in turn, allows to im- 
plement a complete reconstruction of the network state 
by solely measuring the qubit, as shown in Sec. IIIII The 
effects of the major sources of noise will be taken into 
account in Sec. [IV] (see also Appendices E] and [B|. We 
conclude the paper by considering the example of a lin- 
ear chain of oscillators (Sec. |V| and discussing possible 
implementations of our scheme fSec. IVIj) . 



I. HAMILTONIAN AND TIME EVOLUTION 

We consider a network-qubit system whose total 
Hamiltonian at t > is given (in a frame rotating with 
the free Hamiltonian of the qubit) by 

H{t)=Ho + Hi^tit), (1) 

N 

Ho = ^ w„ajja„ -I- ^ Jnm {analn + aiarn) + 

n=l n<m 

+ ^ K„rn {a„ajn + aj^aj^) , (2) 

n<m 

Hintit) = 9{t)a3 (ai + (3) 



where N is the number of oscillators, a„ the bosonic an- 
nihilation operator for the n-th oscillator, a;„ the cor- 
responding local frequency, Jnm and Knm the interac- 
tion strengths between the n-th and the m-th oscillator, 
and g{t) the time- varying coupling strength between the 
qubit and a single oscillator of the network, which we la- 
bel n = 1. The operator is a generic Pauli operator 
for the qubit, belonging to a right-handed tern cti, 1725 cs, 
with [ai,aj] = 2i eijk<Jk- The mutual interactions be- 
tween the oscillators and the qubit in Eqs. ([2]) and ^ 
suggest that information can propagate from any node 
of the network to the first one (and the qubit), in turn 
permitting the reconstruction of the whole network state 
by accessing only the first node. However, these mutual 
interactions also yield the dispersion of any signal along 
the network. Thus, before giving the explicit reconstruc- 
tion method, we first have to solve the inverse problem 
of unravelling the intricate dynamics of this interacting 
system. 




FIG. 1: Graphical representation of the Hamiltonian of 
Eq. Ilj. A qubit (two-level system) is tunably coupled [g{t)] 
to a single constituent (ai) of an oscillator network. The con- 
stituents of the network (aj, with j — 1,...,N) are in turn 
harmonically coupled {Jn,m, Kn^rn)- 



A. Normal modes decomposition 

To study the time evolution of the system, it is con- 
venient to first express the network Hamiltonian in the 
diagonal form 



Hn 



N 



k=l 



(4) 



where h^s are the normal mode operators with corre- 
sponding eigenfrequencies 
monic network is stable, i.e. 
normal modes b = (61,..., 6^) are related to the local 
modes a = (oi, a^v) via [23[ 



We assume that the har- 
Vk > ^ for any k. The 
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a* 



(5) 



where 5 is a 2A^ x 27V symplectic matrix That is, S 
is a transformation that preserves the canonical commu- 
tation relations. It proves convenient to decompose S in 
four N X N blocks. Looking at Eq. ([5]), we see that it is 
possible to write down 



S 



Si S2 



(6) 



where * indicates element-wise complex conjugation (as 
opposed to hermitian conjugation, where the matrix is 
also transposed) . The preservation of bosonic cummuta- 
tion relations imply the constraints 



sisi - {sis2r = 1, 

sls2^{slsiy. 

As a consequence, the inverse of S is given by 



s- 



Si 

-S2 s^ 



(7) 
(8) 



(9) 
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FIG. 2: In the normal modes representation of the oscillator 
network, the qubit is interacting with those modes bk such 
that Gfc 7^ 0. Each mode bk behaves as a simple harmonic os- 
cillator of frequency Uk , and does not interact with the others 
[see Eqs. HO]) and JH}]. If the assumptions (Al) and (A2) 
are verified, the qubit interacts with all the normal modes, 
and can distinguish each mode by its frequency. 

From this, we can express the interaction Hamiltonian of 
Eq. ([3]) in the new basis: 

N 

iJi„t(0 = g{t)a3 (Ckbk + G^fol) > (10) 
fc=i 

Gk = (S, - S,)l,. (11) 

We note that, in the new representation, the qubit in- 
teracts with all the modes bk such that Gk ^ 0. Fig. [2] 
shows a graphical representation of the Hamiltonian ([1]), 
in terms of the normal modes of the oscillator network. 



B. Time evolution for a closed system 

It is convenient to evaluate the time evolutor in the 
normal modes basis. The Hamiltonian ([T]), in an interac- 
tion picture with respect to Hq, can be recasted as 

N 

Hi{t)^Y.^k{tl (12) 
fc=i 

hk{t) = g{t)o^ (Ckhke-''"'* + Glble''"'') . (13) 

We can sec that [hkit),hk'{t')] = if fc 7^ k' . Therefore 
the time evolutor must be of the form Ui{t) = <E}kUk{t), 
where Uk obeys the Schrodingcr equation 

ilk = -ihkUk, (14) 
u(0) = 1. (15) 

We can impose the ansatz Uk = e**'' exp{cr3[/3fe6|, — 
/3^6fe]}, where 4>k and j5k are functions of time. This 
leads to Pk = —w{t)Gl.e^'^''* plus an equation for (j>k that 
we do not need to solve, since the product Jlfc^"^''''*'' is 
just a global phase factor that can be ignored. Then, the 
time evolutor of the system can be given in the closed 
form: 

Uiit) = exp {a3[b^/3(5, t) - f3{g, t)^b] } , (16) 



where (3{g,t) = {I3i{g,t), j3N{g,t)) and 

pk{g,t) = -iG*k ^dsg(s)e^''^^ (17) 
"'0 



II. REALIZING ARBITRARY 
QUBIT-CONTROLLED DISPLACEMENTS 

The time evolutor of Eq. p6p yields a qubit- controlled 
multimode displacement for the bosonic modes bk, char- 
acterized by displacement parameters ±Pk{g,t), the sign 
being determined by the eigenvalue of a^. Notice that 
the displacement (3{g, t) is a functional of the coupling 
strength g{s). The ability to tunc the latter will be cru- 
cial in reconstructing the state of the network. In terms 
of the local modes a, the time evolutor reads (for brevity, 
from now on we will omit the explicit dependence of the 
displacement on g and t) 

C//(i) = exp {0-3(0^0; — cc^a)} , (18) 

where 

(T)=^'(t')- 

where S'^ indicates the transpose of S, while a. = 
(ai, a^)- At this point we make two assumptions: 

• (Al) all the coefficients Gk are different from zero 

• (A2) the normal modes spectrum {vi, ...,vn} is 
non-degenerate 

In physical terms, they imply that the probe qubit in- 
teracts with and can resolve all normal modes bk [see 
Eqs. (fT2TT3l) and Fig. [2]. 

These assumptions are satisfied by generic networks 
(i.e., networks without special symmetries) and, in par- 
ticular, by a linear chain of oscillators. When (Al) and 
(A2) are verified, it becomes possible to assign arbitrary 
values to the displacement vector a, just by controlling 
the length of the interaction time t and by appropriately 
tailoring the time dependence of the coupling g{s). To 
see that this is possible, suppose that we wish to apply 
the operator of Eq. (fT8|) . with generic a of our choice. 
The corresponding vector /3 that has to be applied to the 
normal modes b is given by 

( "J ) - (^')- ( 7' ) ■ (») 

Thanks to the assumption (Al), we can impose an inter- 
action strength profile of the form 
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with Bi coefficients to be determined. Inserting the above 
expression in Eq. ()f 7p . and defining B = (_Bi, B^r), we 
obtain 



-(3* 
/3 



M 



= M 

Ml M2 
M3 Mi 



-B* 
B 



(22) 
(23) 



where M is a 2N x 2N matrix, whose four N x N blocks 
are: 





Gl t Jq 


(24) 






(25) 


{M3)kl 




(26) 


{Mi)M 




(27) 



Matrix M is invertible for long enough interaction times. 
In fact, it is easy to see that assumption (A2) implies 



lim AI = I2N 

i— f 00 



lim det M = 1 (28) 



where I2N is the 2N x 2N identity matrix. This im- 
plies that there must be an interaction time to such that 
detM > for t > to. Inverting Eq. (|22|). we obtain the 
required values of {Bi, Bn): 



-B* 
B 



M- 



-f3* 
/3 



{s''My 



a 
a 



(29) 

A simple Fourier argument as well as the examples we 
studied numerically suggest that the matrix M is invert- 
ible when we take 



t > 



(30) 



that is, when the interaction time is sufficiently long to 
resolve the smallest frequency difference of the system. 
In practice, if we require the components of the displace- 
ment vector OL (or (3) to be large, it might be necessary to 
pick even longer interaction times, since the magnitude 
\g{s)\ is often limited to a maximum value in realistic 
implementations. In fact, one can see from Eqs. (PT|) 
and ((29|) that, keeping a (or /3) fixed, longer interaction 
times imply a smaller amplitude of the coupling strength 
g{s). On the other hand, the interaction time t has to 
be kept small compared to the decoherence timescales of 
the system. The combination of these two facts imposes 
practical limits to the maximum attainable values of |a„| 
(or |/3fc|), thus reducing the extent of the accessible region 
in the phase-space of the oscillator network. 



III. QUANTUM STATE RECONSTRUCTION 

Let us show that the ability to tune the qubit-network 
coupling allows to reconstruct an arbitrary initial state 



of the network just by performing measurements on the 
qubit. To begin, we initialize the system in the state 



Ptot(0) = 



(31) 



where p is the unknown state of the network at t = that 
we want to reconstruct, |-|-) = :^(|<7) + |e)) is the positive 
eigenstate of (Ti, |e) and \g) being respectively the posi- 
tive and negative eigenstates of a^,. We then choose an 
interaction time t > to and a set of local displacement pa- 
rameters a = — ^/2 [with ^ = {£,1, ■■■,(,n) and the factor 
— i included for later convenience], so that a specific pro- 
file of g{t) is determined according to Eqs. (pijl and (^5]) . 
After the interaction, the system has evolved to a state 
Ptot(^) = Ui{t)ptot{Q)Ui{t)^ . We then measure either the 
qubit observable cti or (72, and repeat the experiment a 
sufficient number of times to estimate the average values 
(aj) — tv {ptotit)aj}. By explicit calculation, we get 



with 



{ai)+i{<J2) ^xiO, (32) 



X(6 = tr{pexp($.at-r -a)} (33) 



being the Weyl characteristic function [l5| of the oscil- 
lator network. By repeating the procedure for different 
points ^ of the network phase space, the full character- 
istic function can be measured. We recall that the lat- 
ter gives a complete description of the state of a multi- 
mode systern, equivalent to its Wigner function or den- 
sity matrix [l5|. In contrast with standard tomographic 
reconstructions @, Eq. ([5^ provides a direct link be- 
tween x(€) and the measured data, without the need 
of any integral transform of the latter. In a sense, the 
post-processing typical of quantum tomography is here 
replaced by the pre-processing needed to determine g{s). 
The advantage is that, while the former is performed on 
noisy state-dependent data, the latter involves only the 
state-independent Hamiltonian parameters. As typical 
for any infinite dimensional system, the full reconstruc- 
tion of the state p [i.e., the entire x(^)] is impractical. 
However, a number of interesting properties can be ac- 
cessed given only a finite collection of x(€) values, as we 
illustrate in the following. 



A. Quantum properties without full reconstruction 

The nonclassicality of a continuous variable state is 
generally associated with its Wigner function being neg- 
ative. In turn, a method to probe this nonclassicality 
criteria directly from a finite collection of characteristic 
function values has been recently put forward and exper- 
imentally tested on a single-mode radiation state [16|. 
Our reconstruction scheme is in this respect especially 
suited, providing directly x(^) from measurements. In 
particular, it might open the way to directly estimate 
nonclassicality for multi-mode states of massive oscilla- 
tors. Other nonclassicality criteria, relying on constraints 
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for the P-function, are also testable directly from the 
characteristic function [l^. In addition, also entangle- 
ment can be similarly estimated. In fact, as suggested 
in Ref. [l^, the method outlined there can readily be 
extended to provide lower bounds for entanglement mea- 
sures in the multimode setting. 

Dealing with the characteristic function offers other 
relevant features. For example, it is often the case that 
one is interested in a block of a system (i.e., its re- 
duced states), rather than the whole system — e.g., to 
test entropy-area laws for many-body ground states Q- 
Given x(€)) this can be readily done, since tracing away 
a mode aj simply corresponds to evaluating x(^l^i = 0) 
[see Eq. ([55)) ]. More generally, correlation functions are 
of broad interest — e.g., in many-body model simulations. 
One can apply polynomial or functional (e.g. Gaussian) 
fits to a finite set of characteristic function values, mea- 
sured in the vicinity of the phase-space origin, to esti- 
mate low order moments of the modes a (in dealing with 
noisy data, this approach is preferable to the extraction 
of moments by derivatives [l5|). From those, any corre- 
lation function of the same order can be calculated. A 
particular but relevant case appears when one considers 
Gaussian states. Then, only second moments are neces- 
sary to reconstruct the state and properties thereof 
Also, deviation from Gaussianity can be addressed by 
considering higher order moments. 



B. Temperature measurements 

Let us consider a concrete example of the reconstruc- 
tion method that might be relevant in first experimental 
implementations. Suppose that one expects p to be in a 
thermal state of the Hamiltonian of Eq. ([5]), and wants 
to test this hypothesis. Being thermal states diagonal in 
the normal modes 6, it is convenient to reconstruct the 
characteristic function directly in terms of the latter — 
something that can easily be done using the method out- 
lined above. In the normal mode basis, the characteristic 
function of a generic thermal state is 



IV. NOISE AND ERRORS 

Let us discuss some of the main sources of error that 
could affect our scheme. Firstly, there is the unavoidable 
coupling of the system to the external environment, giv- 
ing rise to decoherence. If this effect can be modelled via 
a standard Markovian master equation, the state of the 
network can still be reconstructed in full detail, at the 
expense of collecting larger amounts of statistical data. 
Secondly, systematic errors might limit the precision to 
which we can control the coupling g(t), meaning that the 
actual displacement parameters will be slightly different 
from the desired values. This will effectively limit the 
phase-space resolution of the reconstructed state. These 
two important sources of error are discussed in detail in 
the sections below. 

Another source of error arises from the experimental 
uncertainties in the Hamiltonian parameters in Eq. ([2|). 
and it affects every stage of our protocol through stan- 
dard error propagation. It is then crucial for the assump- 
tions (Al) and (A2) to be verified for the whole range 
of parameters inside the error bars. If this condition is 
met, the inversion of matrix M in Eq. ([29]) remains well 
defined. Moreover, Eq. ([^5)) guarantees that, for long 
enough interaction times, the uncertainty on will 
be of the same order of the uncertainty on M. 

Finally, as common in many reconstruction protocols, 
errors in the measured data can yield a non-physical re- 
constructed state. In our case, the crucial issue is to 
check whether a finite collection of measured characteris- 
tic function values (with associated uncertainties) is com- 
patible with a positive semidefinite density matrix [nor- 
malization can be satisfied simply by imposing x(0) = !]■ 
This problem can be addressed directly by making use of 
the quantum Bochncr Theorem [isj . for exaraple by us- 
ing the numerical methods developed in Ref. [l6[. There 
it is shown how, by using semidefinite programming, it is 
possible to output a set of characteristic function values, 
compatible with the data but devoid of non-physicalities. 



A. Markovian Decoherence 



N 



Xriv) exp<^ 



(34) 



where A/'(i^A:) = l/(e'^''/-^ — 1) is the number of bosonic ex- 
citations at frequency and temperature T . Following 
the procedure above, one can choose a set of displacement 
parameters /3 = — ry/2 [see Eq. ([^ ]. so that a finite col- 
lection of x(^) values can be reconstructed directly in the 
normal mode basis. Then, standard statistical methods 
can be employed both to test the validity of the thermal 
hypothesis and estimate T. 



To treat environmental noise in our model, it is con- 
venient to work in the normal modes basis of the os- 
cillator network. To simplify matters, we will restrict 
the discussion to the decoherence of the oscillators being 
"diagonal" in terms the normal modes 6^. This can be a 
good model of decoherence when the environmental noise 
is completely uncorrelated between different nodes of the 
network, and when the inter-oscillator couplings Knm are 
small with respect to the local frequencies w„,a;„i (see 
Appendix [Xj for a more detailed discussion) A widely ap- 
plicable model of Markovian decoherence for both the 
qubit and the oscillators is given by the master equation 
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m 

Ptot = -i[Hi{t),ptot] + ^ LkPtot + Sptot- (35) 

k 

The terms responsible for decoherence are 

Lk = ^{Mk + l)V[bk] + ^MkV[bll (36) 

Q = Y(A/', + + YAA,P[a+] + ^P[a3], (37) 

where a+ = |e)(g|,(j_ = |5)(e|, Kk {k = l,...,N) is the 
couphng of each normal mode to the environment (they 
might be in general different, as each normal mode has 
a different frequency), while A/fc = M{i'k) is the thermal 
occupation of the environment at frequency Vk- Fi and 
r2 are the qubit couplings to the environment, the first 
being responsible for thermalization, Mq = N{uiq) being 
the thermal occupation of the environment at frequency 
LOq, while T2 is the strength of additional dephasing mech- 
anisms. Finally, the action of the superopcrator I? on a 
generic operator A is 

V[A]ptot = 2AptotA^ - AUptot - PtotA^A. (38) 

Note that, by using the above model of decoherence, we 
have the simplification that the environment does not 
induce any coupling between the normal modes of the 
network, (see Appendix [A|) . 

By solving the dynamics analytically, it can be shown 
(see Appendix |B]) that Eq. ([5^ has to be modified as 
follows: 

{a{)+z{a2)^x{'n)e~f^^-''\ (39) 
/ being a positive function: 

/(.g, t) = ^t+J2 [Afe (1 - e~^^') \^ik{t)\^ + Tfc(i)] , (40) 

k 

where the explicit forms for p,k and are given by 
Eqs. (|B12p and (|B13p respectively. We can see that the 
effect of decoherence is twofold. Firstly the matrix M, 
which in this case gives (— ^J*,??) = — 2Af(— B*,B), is 
now given by 

(A/i)fci = ^7 ^dse-^(''^-'''>-^^ (41) 

tr/ t Jo 

{M2),k = %\ f dse-'^'"'+'''>-'^' , (42) 
^/ * Jo 

while {M3)ki = {J^'U)ki = (A^i)fe;- A crucial point 

in our protocol is the invertibility of the matrix M, since 
it allows us to assign arbitrary values to r}. We have al- 
ready seen that if the time t is chosen large enough, this 
matrix can be inverted in the C3jS6 Atfc — 0, i.e., detM ^ 0. 
Due to the continuity of the determinant, if the condi- 
tion t <C 1/Kfc, (fc — 1, ...jN) can be verified, then M is 
just slightly perturbed when Kk 7^ 0, so that its determi- 
nant remains different from zero. In cases where t and 



I/kj are of the same order, the invertibility of the matrix 
M should be verified numerically. The second effect of 
decoherence is the appearance of the damping term e^^ 
in Eq. (j39p . meaning that the measured quantity devi- 
ates from the actual value of the characteristic function. 
However, since the function /(g,t) is state-independent 
in our model, and < / < 00, it follows that the right 
hand side of Eq. ([M)) is still a valid representation of the 
quantum state p. This means that we could in principle 
recover the value of %(''?) j if the decoherence parameters 
of the system are known with sufficient accuracy, simply 
by multiplying the measured data by e-^ . Note however 
that this operation also applies to the experimental un- 
certainty, so that an error d in the measured data implies 
a larger error Se-^ in the knowledge of x(*7)- To compen- 
sate for this, the expectation values of the Pauli operators 
in Eq. p9|) have to be measured with higher accuracy 
compared to the decoherence-free case, which necessarily 
requires a larger number of experimental repetitions. Fi- 
nally, we recall that we are considering x in the normal 
modes basis, so the complex vector rj is related to the 
local modes vector ^ via (— C*,C) ~ S'^ {—rj* ,t]). 



B. Noise in the tunable coupling 

Since our protocol relies heavily on the controllability 
of the time-dependent coupling g{t), it is natural to ques- 
tion its robustness against imprecisions in such control. 
We model systematic errors in the controllable coupling 
by substituting 

9{s)^9{s) + as), (43) 

where C is a small white noise component: 

= 0, (44) 
C(ii)C(t2) = e<5(ti-t2), (45) 

and indicates the ensemble average, that is, the aver- 
age over many realizations of the noise [physically, over 
many repetitions of the experiment, where the deter- 
ministic component g{s) is kept fixed]. In Eq. (j45|) . e 
has the dimensions of a coupling strength variance per 
unit of frequency (with our choice of units, this amounts 
to the dimensions of a frequency), and it represents the 
strength of the white noise. Loosely speaking, e repre- 
sents the "thickness" of the curve {s,g{s)). As a result 
of Eq. ((43|) . at each repetition of the experiment the dis- 
placement parameters in the time evolutor of Eqs. (jl6p 
and (fT8|) deviate from the desired values by a small ran- 
dom amount. This yields a finite resolution in our power 
of observation of the oscillator phase-space, meaning that 
we will only be able to measure a coarse-grained version 
of the Characteristic Function, where sub-resolution fea- 
tures are washed out. For simplicity, let us compute this 
phase-space resolution in the normal modes basis. If we 
plug Eq. (231) into Eq. p7)) . we see that the displacement 



parameters f3k in the time evolutor (jl6p are replaced by 
Pk^f3k + Sl3k, (46) 
where S^k are the complex random variables 



<5/3fc = -iGl / dsC(s)e*^'=^ 
Jo 



(47) 



Due to Eqs. (|44|) and (j45|) . we can see that their means 
and covariances are respectively 



0, 



(48) 



(49) 



By diagonalizing V{6l3)^ we can find a new basis for the 
phase-space, such that the noise along different (orthogo- 
nal) directions is uncorrelated. Note how such diagonal- 
ization is independent on the noise strength e, and it only 
depends on the interaction time t and the structure of the 
network. The square roots of the eigenvalues of V{6(3) 
can then be used as a measure of the achievable phase- 
space resolution, along the phase-space directions defined 
by the new basis. When the interaction time is large 
enough, one can see that the diagonal terms [V{Sf3)]kk 
are dominant |24|. Thus, to have an estimate of the order 
of magnitude of our accuracy in phase-space, we can look 
at the diagonal elements of the covariance matrix, which 
have the simple form: 



[V{S(3)] 



kk 



(50) 



It follows that, 0jS 0, first approximation, and assuming 
that the presence of systematic noise in the coupling 
strength is well modelled by Eq. (|43p . the smallest resolv- 
able feature in the phase space of the oscillator network 
scales only sublinearly with e and t: 



\6/3k 



\Gk\vrt. 



(51) 



V. EXAMPLE: LINEAR CHAIN WITH 
CONSTANT COUPLINGS 

Let us give a concrete example of a quantum oscillator 
network where the presented ideas can be applied. Con- 
sider a linear chain of N oscillators, each having the same 
local frequency w, and where only the nearest-neighbours 
interact with a coupling strength J„,„+i = Kn,n+i = </, 
constant along the chain. The J's are often referred to as 
hoppings. We assume that the qubit is tunably coupled 
to the first oscillator of the chain. The system is sketched 
in Fig. [21 In terms of the parameters appearing in the 
Hamiltonian Hq [see Eq. ([2])], we have 



I^nm — Jnrn — 







n<N, 
otherwise. 



(52) 
(53) 



ai 



02 
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FIG. 3: Sketch of a qubit tunably coupled to a linear chain 
of oscillators with constant nearest-neighbour couplings. 



The diagonalization of the Hamiltonian Hq can be per- 
formed analytically, yielding the spectrum (here k = 
1,...,7V) 



Uk = \/oj{uj + 2efe), 
, Trfc 

2 J cos 



iV + 1 



(54) 
(55) 



One can see that the above spectrum is non-degenerate, 
thanks to the fact that the cosine is monotone in the in- 
terval [0,7r]. Thus, the linear chain of oscillators verifies 
the assumption (A2). The symplectic matrix 5, connect- 
ing the local modes a„ to the normal modes bk, can also 
be expressed in analytical form. Its main blocks 5*1 and 
S2 are (see Section [O]) 



iSl)kn - 
{S2)kn = - 

Tfc = tanh" 



iV+ 1 



cosh (rfc) sin 



iV+ 1 



sinh (rfc) sin 



T:kn 
N + 1 
Trkn 



£k 



- £k + . 

If we combine Eqs. (f5S)l . ([57)) and pT|) . we have 
^ / " ^. . I '^k 



iV+ 1 



N + 1 



(56) 
(57) 
(58) 

(59) 



which is different from zero for any /c € (1, N). There- 
fore, also assumption (Al) is verified. Thus, the quan- 
tum state of a linear chain of oscillators with constant 
nearest-neighbour couplings can be fully reconstructed, 
by using a single qubit coupled to one end of the chain. 
As an example. Fig. [4| shows some quantities of interest 
for the reconstruction protocol of a linear chain of = 8 
oscillators. 



VI. POSSIBLE EXPERIMENTAL 
IMPLEMENTATIONS 

The reconstruction scheme described here is based on 
a rather ubiquitous dynamics. Essentially, it requires a 
harmonic coupling between the oscillators and a bilinear 
qubit-oscillator interaction. An experimental platform 
that is particularly mature for our purposes is given by 
a chain of ions in a linear trap. There the harmonic 
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FIG. 4: Reconstruction protocol for a linear chain of A'^ = 8 
oscillators. In all plots we consider nearest-neighbour cou- 
plings Jnm = Knm = J6m,n+i, with J = 0.2tj, decoher- 
ence rates Kk = k = W~^uj temperature T = 200aj and a 
coupling noise strength e = 10~''aj. Note that we are re- 
quiring a motional quality factor Q ~ ~ 5 x 10^. We 
assumed the decoherence of the qubit to be negligible. Plot 
(a) shows the determinant of the matrix M, which we require 
to be different from zero in our protocol. We see that M 
becomes invertible for t > 50/uj. In plots (b) and (c), we 
have considered a phase-space region |r7fe| < 2 in the normal 
modes basis, and the total interaction time has been fixed to 
t = 100 X 27r/cj. Plot (b) shows a portion of the interaction 
strength profile required to obtain the displacement param- 
eters /3 = (— 1, — 1, — 1), which allows us to reconstruct 
the phase space point t/ = (2, 2,..., 2), lying at the bound- 
ary of the considered region. With our choice of parame- 
ters, it is necessary to access maximal coupling strengths of 
Iplmax — O.OSoj to reach this phase space point. In plot (c), the 
effect of decoherence on the measured value of the character- 
istic function is shown. We fixed ri2 = rjs = ... — rjN = 2 and 
Im{r;i} — 0. The plot shows the quantity as a function 
of Re{r;i}. We see that near the boundary of the considered 
phase-space region the quantity measured via the qubit ex- 
pectation values corresponds to about 20-23% of the actual 
value of the characteristic function [see Eq. (|39[) ]. Plot (d) 
shows the maximum achievable phase-space resolution along 
the directions that diagonalize V{Sf3) [we take the square 
roots of its eigenvalues Xk as a measure of the phase-space 
accuracy]. The asymptotical behaviour for large t is given by 
x/AT — |Gfe|\/rf, as expected. As a final remark, we note that 
the considered range of interaction times is consistent with 
Eq. (|XT8)) . which gives a necessary condition for the valid- 
ity of the master equation used. Indeed, from Eq. (|A18|l one 
can estimate that our treatment is valid for t <^ tmax, where 
W = miufc [^^] f ~ 2 X W^{2tv/lj). 

dynamic is provided by the Coulomb interaction, as re- 
cently demonstrated in Ref. [l^. In addition, the re- 
quired qubit-oscillator coupling is standard [lo| : two elec- 
tronic levels of one ion provide the qubit, whereas the 
coupling is realized via a standing laser wave [l9| . For 
example, consider a linear chain of A'^ = 8 ions, with 
Wn uj, Jn^m = Knrn ^ ^m,ri-f i0.2(.j and take an in- 
teraction time t ^ 100(27r/aj). Then |(7(s)|max ^ O.OSoj 
allows to reconstruct states with x(^) having support in 



|^„| < 2. Notice that a modest motional quality factor 
Q >5x 10"^ is required (see Fig.|4]), and it is sufficient to 
vary the laser power on a time scale of the eigenfrequen- 
cies (typically of the order of MHz) in order to realize 
the desired profile of g(s). We stress again that only one 
ion needs to be illuminated in order to reconstruct the 
motional state of the entire chain. 

Similar couplings can be envisaged also in a circuit 
QED setup [201, where many stripline waveguides can 
be capacitively coupled mimicking a harmonic network, 
and a single superconducting qubit can be coupled to one 
of the resonators. In this setting, the realization of the 
Hamiltonian (jl2p requires a regime where the qubit level 
split ting is negligible compared to the dipolar coupling. 
Ref. [2l| is promising in this direction, as it demonstrates 
independent tunability of these two parameters. 

Being our model quite general, its implementation can 
be envisaged also in other experimental platforms, such 
as nanomechanical oscillators or microcavities. In fact, 
the recent effort to build complex quantum networks led 
to impressive experimental progresses. In this context, 
the question of extracting information from a quantum 
network by only accessing a limited portion of it is of 
interest from a general viewpoint. We provided here a 
solution to this problem for a network of quantum oscil- 
lators and we expect that further investigations for the 
case of different constituents will be relevant in the fu- 
ture. 
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Appendix A: Derivation of the master equation for 
the oscillator network 

Let us start by writing down a Hamiltonian that in- 
cludes the harmonic oscillator network, its environment, 
and their mutual interaction. We assume that each oscil- 
lator mode a„ is interacting with a local environment, de- 
scribed by a continuum of bosonic modes Cn(aj), and that 
their interaction is bilinear, with frequency-dependent 
coupling strength /n(a;). Thus we take 

^^tot ~ Ho+ Henv + HsE, (Al) 
TV 

n— 1 n<rn 

+ ^ Knm {andrn + olal^) , (A2) 

n<7n 

HsE = + 4) / dLjULo)[cn{io) + cliio)], (A3) 

Henv = ^ / dw UJcl{uj)Cn{i^), (A4) 
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We now assume that the environment is completely un- 
correlated between different points of the oscillator net- 
work, and that the corresponding bosonic modes are in- 
dependent, that is: 



(A5) 



If in Eq. (|Aip we switch to the normal modes bk, which 
diagonalize Hq, we can rewrite the interaction Hamilto- 
nian of Eq. (|A3P as 



HsE = Y.(Kkbk + KA) J dooUiu)[cn{Lo) + ct (c.)], 

kn 



^^{Si- 82)* 



(A6) 
(A7) 



where Si, S2 are the N x N matrices defined in Eq. 
Assuming that the coupling between system and environ- 
ment is weak, we can neglect the counter-rotating terms 

M 

HsE^Y^j ^^Uu:)[Kkhlcn{uj) + ^nkbkC^n{^)]. (A8) 

kn 

We take the further approximation that the environmen- 
tal coupling is the same at any location of the network: 
/„(w) = /(w), so that we can perform the sum over the 
index n, and write 

HsE^Y.! d^/(^)[&fc4(^) + bldk{oJ% (A9) 

k 

dkiuj) ^Y.'^nkCni^)- (AlO) 



Then, the free Hamiltonian of the environment can be 
rewritten as 

i/cnv = YjJ'^^ ujdl{uj)dk{uj) + O . (A15) 

It follows that, at the zeroth order in Knml^qi each 
/c-subspace evolves independently, according to a total 
Hamiltonian 

Ht ~ Vkblbk + j Aujf{u)[budl{io) + hldk{i^)] + 

+ dij Ljdl{oj)dk{uj), (A16) 



where the operators can be treated as bosonic. At 
this point, one can apply standard techniques to derive 
separately the master equation for each normal mode bk 
[15j . Putting together all the modes, one can easily derive 
the oscillator part of the master equation of Eq. ([551) . 
particular, the coupling parameters Kk are given by 



(A17) 



Due to the approximations used to derive it, our master 
equation is only valid for timescales such that 



t <^ min 

k 



1 



kJ 



(AI8) 



where Afk is the number of thermal bosonic excitations 
at frequency Vk- 

Appendix B: Solving the master equation 



Note that the operators dk are not bosonic, since the 
matrix -d is not unitary in general. Indeed, Eq. (jATj 
together with Eqs. ([7| and ^ imply 



t?^t9 = 1 - (SlSi)* - {SIS2)* + '2{SlS2)* 



(All) 



However, if the "active" terms in the Hamiltonian are 
weak, that is |/'Sr„„i| ^ ujn,uj„i [for any combination of 
n,m — see Eq. (|A2p ]. then one can see that 



S2 



(A12) 



where the expression O (■^) indicates the order of mag- 
nitude of the ratios K^m/'^n, K^m/'^m- It follows that 
1? is approximately unitary and the operators dk become 
approximately bosonic. [25[ 



K 



[dk{uj),dl,{u:')]^5kk'5{Lo-uj')+o(^-^^ 



(A13) 



(A14) 



To solve the master equation, we consider a representa- 
tion in which a matrix of characteristic functions is used 
to describe the state of the coupled system. We decom- 
pose the total density matrix at time t as 

PtotW = Pe{t) ® |e)(e| + pg{t) ® \9){g\ + 

+ p+{t)®\e){g\+p-{t)®\g){el (Bl) 

where the operators pj {j = e,g,+,—) belong to the 
oscillators. If we define the characteristic function for 
each element as 

Xj(/3,i) =tr6,,...,h« {pj{t)cxp{b^(3-f3^b)}, (B2) 

we can write 

x(/3, t) = Xe(/3, m (e| + Xg{f3, t)\g) {g\ + 

+ X+{f3,m{g\+X-{f3,t)\g){e\. (B3) 

In this formalism, the expectation values of the Pauli 
operators in Eq. (15) of the main text yield: 

tr{ptot(i)((Ti + ztT2)} = 2x-{0,t) (B4) 
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It follows that we only need to compute the time evo- 
lution of the element x-- By using standard techniques 
[l5l | , we can convert the bosonic operators bk and b], into 
differential operators for the characteristic function. Ap- 
plying this procedure to our master equation we find that 
X- is decoupled from the other elements, and obeys the 
following equation: 

k 

+ ^£,X--7X-, (B5) 
fe 

Ck = - Y {Pkdh + Pldpi + 2Afe|/3,|2) , (B6) 
Afc = AAfc + i 7 = Ti{Nq + 1/2) + 2T2. (B7) 
The corresponding solution is 

X-{fi,t)=x- (/9(0 + r7(t),0) e-^*x 





k 


m-- 




vit)-- 




Vk{t) 


= 2iGl f\ 
Jo 


tik{t) 


sinh 


Tkit) 


Jo 



dsg{s)e'''^''^ sinh ^s, 



(B8) 

(B9) 
(BIO) 

(Bll) 
(B12) 

(B13) 



The initial state of Eq. ((3T|) of the main text implies 
the initial condition x-(/3iO) = l/2x(/3), where x(/3) 
is the characteristic function of the initial state of the 
oscillators, expressed in the normal modes basis. Then, 
Eq. ((391) follows by substituting /3 = in Eq. (|B8l) 
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